TD6 Traitement de Signal¶

EXERCICE I¶

Q1.1¶

In [13]:
import numpy as np
import matplotlib.pyplot as plt
In [14]:
Fe = 4  #fréquence d'échantillonnage (Hz)
Te = 1 / Fe  #période d'échantillonnage (s)
duration = 10  #durée (s)

#vecteur temporelle
t = np.arange(0, duration, Te)

#definition du signal s(t)
s_t = 3 + np.sin(2 * np.pi * t)

plt.figure(figsize=(10, 4))
plt.plot(t, s_t, label="s(t) = 3 + sin(2πt)")
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal s(t) sur 10 secondes")
plt.legend()
plt.grid()
plt.show()
No description has been provided for this image

Q1.2¶

In [15]:
from scipy.signal import lfilter
In [16]:
help(lfilter)
Help on function lfilter in module scipy.signal._signaltools:

lfilter(b, a, x, axis=-1, zi=None)
    Filter data along one-dimension with an IIR or FIR filter.

    Filter a data sequence, `x`, using a digital filter.  This works for many
    fundamental data types (including Object type).  The filter is a direct
    form II transposed implementation of the standard difference equation
    (see Notes).

    The function `sosfilt` (and filter design using ``output='sos'``) should be
    preferred over `lfilter` for most filtering tasks, as second-order sections
    have fewer numerical problems.

    Parameters
    ----------
    b : array_like
        The numerator coefficient vector in a 1-D sequence.
    a : array_like
        The denominator coefficient vector in a 1-D sequence.  If ``a[0]``
        is not 1, then both `a` and `b` are normalized by ``a[0]``.
    x : array_like
        An N-dimensional input array.
    axis : int, optional
        The axis of the input data array along which to apply the
        linear filter. The filter is applied to each subarray along
        this axis.  Default is -1.
    zi : array_like, optional
        Initial conditions for the filter delays.  It is a vector
        (or array of vectors for an N-dimensional input) of length
        ``max(len(a), len(b)) - 1``.  If `zi` is None or is not given then
        initial rest is assumed.  See `lfiltic` for more information.

    Returns
    -------
    y : array
        The output of the digital filter.
    zf : array, optional
        If `zi` is None, this is not returned, otherwise, `zf` holds the
        final filter delay values.

    See Also
    --------
    lfiltic : Construct initial conditions for `lfilter`.
    lfilter_zi : Compute initial state (steady state of step response) for
                 `lfilter`.
    filtfilt : A forward-backward filter, to obtain a filter with zero phase.
    savgol_filter : A Savitzky-Golay filter.
    sosfilt: Filter data using cascaded second-order sections.
    sosfiltfilt: A forward-backward filter using second-order sections.

    Notes
    -----
    The filter function is implemented as a direct II transposed structure.
    This means that the filter implements::

       a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[M]*x[n-M]
                             - a[1]*y[n-1] - ... - a[N]*y[n-N]

    where `M` is the degree of the numerator, `N` is the degree of the
    denominator, and `n` is the sample number.  It is implemented using
    the following difference equations (assuming M = N)::

         a[0]*y[n] = b[0] * x[n]               + d[0][n-1]
           d[0][n] = b[1] * x[n] - a[1] * y[n] + d[1][n-1]
           d[1][n] = b[2] * x[n] - a[2] * y[n] + d[2][n-1]
         ...
         d[N-2][n] = b[N-1]*x[n] - a[N-1]*y[n] + d[N-1][n-1]
         d[N-1][n] = b[N] * x[n] - a[N] * y[n]

    where `d` are the state variables.

    The rational transfer function describing this filter in the
    z-transform domain is::

                             -1              -M
                 b[0] + b[1]z  + ... + b[M] z
         Y(z) = -------------------------------- X(z)
                             -1              -N
                 a[0] + a[1]z  + ... + a[N] z

    Examples
    --------
    Generate a noisy signal to be filtered:

    >>> import numpy as np
    >>> from scipy import signal
    >>> import matplotlib.pyplot as plt
    >>> rng = np.random.default_rng()
    >>> t = np.linspace(-1, 1, 201)
    >>> x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1) +
    ...      0.1*np.sin(2*np.pi*1.25*t + 1) +
    ...      0.18*np.cos(2*np.pi*3.85*t))
    >>> xn = x + rng.standard_normal(len(t)) * 0.08

    Create an order 3 lowpass butterworth filter:

    >>> b, a = signal.butter(3, 0.05)

    Apply the filter to xn.  Use lfilter_zi to choose the initial condition of
    the filter:

    >>> zi = signal.lfilter_zi(b, a)
    >>> z, _ = signal.lfilter(b, a, xn, zi=zi*xn[0])

    Apply the filter again, to have a result filtered at an order the same as
    filtfilt:

    >>> z2, _ = signal.lfilter(b, a, z, zi=zi*z[0])

    Use filtfilt to apply the filter:

    >>> y = signal.filtfilt(b, a, xn)

    Plot the original signal and the various filtered versions:

    >>> plt.figure
    >>> plt.plot(t, xn, 'b', alpha=0.75)
    >>> plt.plot(t, z, 'r--', t, z2, 'r', t, y, 'k')
    >>> plt.legend(('noisy signal', 'lfilter, once', 'lfilter, twice',
    ...             'filtfilt'), loc='best')
    >>> plt.grid(True)
    >>> plt.show()

Q1.3¶

In [21]:
#coeffs du filtre moyenneur (fenêtre de 4 points)
coeffs_moyenneur = np.ones(4) / 4

#application du filtre moyenneur sur signal s(t)
s_moyenne = lfilter(coeffs_moyenneur, 1, s_t)

plt.figure(figsize=(10, 4))
plt.plot(t, s_t, label="Signal original s(t)")
plt.plot(t, s_moyenne, label="Signal filtré (moyenneur 4 points)", linestyle='--', color='r')
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Filtre moyenneur 4 points")
plt.legend()
plt.grid()
plt.show()
No description has been provided for this image

Q1.4¶

In [9]:
#def du filtre dérivateur : différencier le signal par différence finie
coeffs_derivateur = [1, -1]

#application du filtre dérivateur au signal s(t)
s_deriv = lfilter(coeffs_derivateur, 1, s_t) / Te  #divise par Te pour obtenir l'approximation de la dérivée

plt.figure(figsize=(10, 4))
plt.plot(t, s_t, label="Signal original s(t)")
plt.plot(t, s_deriv, label="Signal dérivé (filtre dérivateur)", linestyle='--')
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Filtre dérivateur")
plt.legend()
plt.grid()
plt.show()
No description has been provided for this image

Q1.5¶

In [20]:
#calcul de la dérivée théorique
#dérive théorique de 3*sin(2pi*t) est s'(t) = 2pi*cos(2pi*t)
s_deriv_theorique = 2 * np.pi * np.cos(2 * np.pi * t)

plt.figure(figsize=(10, 4))
plt.plot(t, s_t, label="Signal original s(t)")
plt.plot(t, s_deriv, label="Dérivée par filtre dérivateur", linestyle='--')
plt.plot(t, s_deriv_theorique, label="Dérivée théorique s'(t) = 2πcos(2πt)")
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Comparaison entre la dérivée théorique et la dérivée par filtre")
plt.legend()
plt.grid()
plt.show()
No description has been provided for this image

Conclusion¶

La dérivée obtenue par le filtre dérivateur est une approximation de la dérivée théorique, mais elle présente des différences dues à la méthode de différence finie utilisée. Le filtre dérivateur est un bon moyen pour estimer la dérivée, mais il peut être affecté par des décalages et des erreurs liées à la précision de l'approximation. Notamment ici avec un filtre dérivateur simple.

EXERCICE II¶

Q2.1¶

In [25]:
from scipy.io import wavfile
In [26]:
import plotly.graph_objects as go
In [23]:
#load file
sample_rate, signal = wavfile.read("Note_Mi.wav")
#normaliser si le signal est en entier
signal = signal / np.max(np.abs(signal))

plt.figure(figsize=(10, 4))
plt.plot(np.arange(len(signal)) / sample_rate, signal)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal de la note Mi")
plt.grid()
plt.show()
/var/folders/2l/zssv4k554m3gbdphlp5n6dtm0000gn/T/ipykernel_1073/1160791914.py:2: WavFileWarning: Chunk (non-data) not understood, skipping it.
  sample_rate, signal = wavfile.read("Note_Mi.wav")
No description has been provided for this image
In [27]:
fig = go.Figure()

#ajout du signal audio à la figure
fig.add_trace(go.Scatter(x=np.arange(len(signal)) / sample_rate, y=signal, mode='lines', name='Signal Audio'))

#config graph
fig.update_layout(
    title='Signal Audio Interactif',
    xaxis_title='Temps (s)',
    yaxis_title='Amplitude',
    hovermode='x',  #affiche valeurs en fonction de la position de la souris
)

fig.show()

Q2.2¶

In [28]:
from scipy.fft import fft
In [30]:
#calcul de la TFD et module
signal_fft = fft(signal)
N = len(signal)
frequencies = np.fft.fftfreq(N, 1/sample_rate)
magnitude = np.abs(signal_fft) / N  #module de la TFD normalisé
phase = np.angle(signal_fft) #argument de la TFD normalis
In [31]:
print(signal_fft)
[-12.73790814-0.j          -5.58583009+2.49197276j
   0.88471877-2.93342473j ...  -9.95848843+0.5295149j
   0.88471877+2.93342473j  -5.58583009-2.49197276j]

Q2.3¶

In [44]:
#spectre monolatéral
positive_freqs = frequencies[:N//2]
positive_magnitude = magnitude[:N//2]
positive_phase = phase[:N // 2]

plt.figure(figsize=(12, 6))

#spectre de magnitude
plt.subplot(2, 1, 1)
plt.plot(positive_freqs, positive_magnitude)
plt.xlabel("Fréquence (Hz)")
plt.ylabel("Amplitude")
plt.title("Spectre de magnitude (monolatéral) du signal")
plt.grid()

#spectre de phase
plt.subplot(2, 1, 2)
plt.plot(positive_freqs, positive_phase, color='orange')
plt.xlabel("Fréquence (Hz)")
plt.ylabel("Phase (radians)")
plt.title("Spectre de phase (monolatéral) du signal")
plt.grid()

plt.tight_layout()
plt.show()
No description has been provided for this image

Q2.4¶

In [34]:
#ajouter bruit blanc
noise_amplitude = 0.2  #ajustez selon le niveau de bruit désiré
noisy_signal = signal + noise_amplitude * np.random.normal(size=len(signal))

plt.figure(figsize=(10, 4))
plt.plot(np.arange(len(noisy_signal)) / sample_rate, noisy_signal)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal bruité (note Mi avec bruit blanc)")
plt.grid()
plt.show()
No description has been provided for this image
In [ ]:
 

Q2.5¶

In [35]:
from scipy.signal import butter, lfilter
In [46]:
#paramètres du filtre Butterworth (ordre et fréquence de coupure)
ordre_iir = 5  #ordre du filtre
fc = 150  #fréquence de coupure en Hz (ajustez en fonction du spectre)

#conception du filtre passe-bas
b, a = butter(ordre_iir, fc / (sample_rate / 2), btype='low')

#filtrage du signal bruité
filtered_signal_iir = lfilter(b, a, noisy_signal)

plt.figure(figsize=(10, 4))
plt.plot(np.arange(len(filtered_signal_iir)) / sample_rate, filtered_signal_iir)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal filtré avec un filtre IIR")
plt.grid()
plt.show()
No description has been provided for this image
In [ ]:
 

Q2.6¶

Pour choisir un ordre et une fréquence de coupure optimale pour le filtre, examinons d'abord les caractéristiques du signal pour déterminer :

  1. Identifier la fréquence fondamentale de la note "Mi" Ici c'est un Mi2 donc on a vu au TD3 -> environ 83Hz (exactement 82.41Hz) La fréquence de la note "Mi" dépend de l'octave : Mi3 est à environ 165 Hz (par exemple), En fonction du contenu du signal et de l'analyse spectrale en question 3, identifions la fréquence dominante dans le spectre de magnitude, qui correspond à la fréquence fondamentale de la note.

  2. Choisir la fréquence de coupure Pour supprimer le bruit sans affecter la note elle-même, il est optimal de choisir une fréquence de coupure légèrement supérieure à la fréquence fondamentale. Cela permet de conserver la fondamentale ainsi que quelques premières harmoniques, souvent souhaitables pour la clarté et la fidélité du son. Ici, 100 - 150 Hz

  3. Choisir l'ordre du filtre Le choix de l'ordre dépend du compromis entre la pente de coupure et la stabilité : Pour un filtre IIR de type Butterworth, un ordre de 4 à 6 est généralement adéquat pour des applications audio. Un ordre plus élevé donnera une transition plus abrupte mais peut introduire des distorsions de phase. Dans la pratique, un ordre de 5 ou 6 fonctionne bien pour ce type de signal et permet une atténuation efficace du bruit tout en préservant les composants de basse fréquence du signal original.

Q2.7¶

In [37]:
from scipy.signal import firwin, convolve
In [38]:
#paramètres du filtre FIR
ordre_fir = 50  #ajuster pour un résultat similaire au filtre IIR
fc_fir = fc  #utiliser la même fréquence de coupure que pour l'IIR

#conception du filtre FIR passe-bas
fir_coeffs = firwin(ordre_fir + 1, fc_fir / (sample_rate / 2))

#application du filtre FIR
filtered_signal_fir = convolve(noisy_signal, fir_coeffs, mode='same')

plt.figure(figsize=(10, 4))
plt.plot(np.arange(len(filtered_signal_fir)) / sample_rate, filtered_signal_fir)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal filtré avec un filtre FIR")
plt.grid()
plt.show()
No description has been provided for this image
In [ ]:
 

Q2.8¶

Un filtre FIR de type passe-bas avec un ordre de 50 et une fréquence de coupure de 1000 Hz est souvent efficace et peut donner une performance similaire au filtre IIR pour des signaux relativement simples comme celui-ci.

Q2.9¶

Les filtres IIR et FIR peuvent tous deux atténuer le bruit. Le filtre IIR est souvent plus efficace pour des calculs rapides, mais peut entraîner des distorsions de phase. Le filtre FIR, quant à lui, assure une réponse en phase linéaire et une meilleure stabilité, au prix de plus de calculs (car un ordre plus élevé est nécessaire pour obtenir un résultat similaire). Dans cette situation, le filtre FIR est préférable pour un signal musical où une réponse en phase linéaire est essentielle.

In [ ]:
 
In [ ]: